In [511]:
import pandas as pd
df = pd.read_csv('data/PCA_Task.csv', header=0, index_col=0)
df.columns = ['SEM', 'EXE', 'GEN']
import numpy as np
import matplotlib.pyplot as plt
import matplotlib.patches as mpatches
from mpl_toolkits.axes_grid1 import make_axes_locatable
import matplotlib
matplotlib.rcParams.update({'font.size': 14})
from hierHeatmap import *
df = hier_clust(df)
fig = plt.figure()
ax = fig.add_subplot(111)
im = ax.matshow(df, vmin=-1, vmax=1, interpolation='nearest', cmap=plt.cm.RdBu_r, aspect='equal')
divider = make_axes_locatable(plt.gca())
cax = divider.append_axes("right", "15%", pad="20%")
plt.colorbar(im, cax=cax)
ax.set_yticks(range(df.shape[0]))
ax.set_xticks(range(df.shape[1]))
ax.xaxis.set_ticks_position('top')
ax.set_xticklabels(list(df.columns), rotation='75', ha='left')
ax.set_yticklabels(list(df.index), size='large')
plt.savefig('Results/Fig2a_PC_tasks.png', dpi=300, transparent=True, bbox_inches='tight')
plt.show()
In [5]:
import numpy as np
import matplotlib.pyplot as plt
import matplotlib
matplotlib.rcParams.update({'font.size': 14})
# Make a random array and then make it positive-definite
eigvals = pd.read_csv('data/PCA_Task_eigen.csv', header=0, index_col=0).values
fig = plt.figure()
sing_vals = np.arange(len(eigvals)) + 1
plt.plot(sing_vals, eigvals, 'ro-', linewidth=2)
plt.plot([3.5, 3.5], [0, 6], ls='--', c='.3')
plt.title('Scree Plot: Phenotypical Variation', size='xx-large')
plt.xlabel('Principal Component', size='x-large')
plt.ylabel('Eigenvalue', size='x-large')
leg = plt.legend(['Eigenvalues from PCA'], loc='best', borderpad=0.3,
shadow=False, prop=matplotlib.font_manager.FontProperties(size='large'),
markerscale=0.4)
leg.get_frame().set_alpha(0.4)
leg.draggable(state=True)
plt.savefig('Results/Fig2a_PC_Eigen.png', dpi=300, transparent=True, bbox_inches='tight')
plt.show()
In [6]:
import pandas as pd
df = pd.read_csv('./data/PCA_QD.csv', header=0, index_col=0)
df.columns = ['AD', 'SOC', 'DYSL', 'ATT']
import numpy as np
import matplotlib.pyplot as plt
from hierHeatmap import *
df = hier_clust(df)
fig = plt.figure()
ax = fig.add_subplot(111)
im = ax.matshow(df, vmin=-1, vmax=1, interpolation='nearest', cmap=plt.cm.RdBu_r, aspect='equal')
divider = make_axes_locatable(plt.gca())
cax = divider.append_axes("right", "15%", pad="20%")
plt.colorbar(im, cax=cax)
ax.set_yticks(range(df.shape[0]))
ax.set_xticks(range(df.shape[1]))
ax.xaxis.set_ticks_position('top')
ax.set_xticklabels(list(df.columns), rotation='75', ha='left')
ax.set_yticklabels(list(df.index), size='large')
plt.savefig('Results/Fig2b_PC_QD.png', dpi=300, transparent=True, bbox_inches='tight')
plt.show()
In [21]:
# Make a random array and then make it positive-definite
eigvals = pd.read_csv('data/PCA_QD_eigen.csv', header=0, index_col=0).values
fig = plt.figure()
sing_vals = np.arange(len(eigvals)) + 1
plt.plot(sing_vals, eigvals, 'ro-', linewidth=2)
plt.plot([4.5, 4.5], [0, 6], ls='--', c='.3')
plt.title('Scree Plot: Questionnaire Items', size='xx-large')
plt.xlabel('Principal Component', size='x-large')
plt.ylabel('Eigenvalue', size='x-large')
#I don't like the default legend so I typically make mine like below, e.g.
#with smaller fonts and a bit transparent so I do not cover up data, and make
#it moveable by the viewer in case upper-right is a bad place for it
leg = plt.legend(['Eigenvalues from PCA'], loc='best', borderpad=0.3,
shadow=False, prop=matplotlib.font_manager.FontProperties(size='large'),
markerscale=0.4)
leg.get_frame().set_alpha(0.4)
leg.draggable(state=True)
plt.savefig('Results/Fig2b_PC_Eigen.png', dpi=300, transparent=True, bbox_inches='tight')
plt.show()
In [3]:
import glob, os
path = glob.glob('.\Results\DMN_16_REVISION\*.npy')
print path
In [4]:
import numpy as np
from os.path import expanduser
region_labels_fn = 'data/data_cross_corr_Bzdok_DMN16_preprocessed_ROIS.npy'
beh_keysfn = 'data/data_MWQ_keys.npy'
keys = np.load(expanduser(beh_keysfn))[1:]
MWQ_keys = np.append(keys, keys)
y_loadings_1 = np.stack((-np.load(path[0])[1][:,0], np.load(path[-1])[1][:,5], np.load(path[1])[1][:,2]),axis=1)
y_loadings_2 = np.stack((np.load(path[0])[1][:,3], np.load(path[-1])[1][:,1], -np.load(path[1])[1][:,3]),axis=1)
x_loadings_1 = np.stack((-np.load(path[0])[0][:,0], np.load(path[-1])[0][:,5], np.load(path[1])[0][:,2]),axis=1)
x_loadings_2 = np.stack((np.load(path[0])[0][:,3], np.load(path[-1])[0][:,1], -np.load(path[1])[0][:,3]),axis=1)
# Transform the flatten weights for network strength back to matrices for visualisation
idx = np.triu_indices(16, 1)
corr_mat_1 = np.zeros((16, 16, 3))
for i in range(3):
this_mat = np.zeros((16, 16))
this_mat[idx] = x_loadings_1[:, i]
corr_mat_1[..., i] = this_mat + this_mat.T
corr_mat_2 = np.zeros((16, 16, 3))
for i in range(3):
this_mat = np.zeros((16, 16))
this_mat[idx] = x_loadings_2[:, i]
corr_mat_2[..., i] = this_mat + this_mat.T
corr_mat_1_flat = np.hstack((corr_mat_1[..., 0], corr_mat_1[..., 1], corr_mat_1[..., 2]))
corr_mat_2_flat = np.hstack((corr_mat_2[..., 0], corr_mat_2[..., 1], corr_mat_2[..., 2]))
x_loadings = np.vstack((corr_mat_1_flat, corr_mat_2_flat))
y_loadings = np.vstack((y_loadings_1, y_loadings_2))
In [5]:
print ' RTS BOOTS'
print y_loadings_1[:, 1:]
In [6]:
n_components = 2
n_connects = 120
n_areas = 16
# slide 2: SCCA summary plots
RS_keys = ['1', '2', '3', '4', '1', '2', '3', '4', '1', '2', '1', '2', '3', '4', '5', '6']
from numpy.ma import masked_array
import matplotlib
import matplotlib.pyplot as plt
from mpl_toolkits.axes_grid1 import make_axes_locatable
mask_brain = np.hstack((np.ones((32,32)), np.zeros((32,16))))
mask_brain_rev = np.array(mask_brain == 0, dtype=np.int)
mask_thoughts = np.hstack((np.ones((26,2)),np.zeros((26,1))))
mask_thoughts_rev = np.array(mask_thoughts == 0, dtype=np.int)
def RSC_plot(mat, ax):
mat_b = masked_array(mat, mask=mask_brain)
mat_a = masked_array(mat, mask=mask_brain_rev)
im_a = ax.matshow(mat_a, vmin=-0.9, vmax=0.9, cmap=plt.cm.RdBu_r)
im_b = ax.matshow(mat_b, vmin=0, vmax=1, cmap=plt.cm.Greys)
ax.locator_params(nbins=3)
ax.set_xticks(np.arange(n_areas*3))
ax.set_xticklabels(RS_keys * 3, fontsize='large')
ax.set_yticks(np.arange(n_areas*2))
ax.set_yticklabels(RS_keys * 2, fontsize='large')
ax.plot([-0.5, 31.5], [-0.5, 31.5], ls='--', c='.3')
ax.plot([15.5,47.5], [-0.5,31.5], ls='--', c='.3')
ax.plot([-0.5,15.5], [15.5,31.5], ls='--', c='.3')
ax.plot([31.5, 47.5], [-0.5,15.5], ls='--', c='.3')
ax.vlines(15.5, -0.5, 31.5)
ax.vlines(31.5, -0.5, 31.5)
ax.hlines(15.5, -0.5, 47.5)
def MWQ_plot(mat, ax):
mat_b = masked_array(mat, mask=mask_thoughts)
mat_a = masked_array(mat, mask=mask_thoughts_rev)
im_a = ax.matshow(mat_a, vmin=-0.9, vmax=0.9, cmap=plt.cm.RdBu_r)
im_b = ax.matshow(mat_b, vmin=0, vmax=1, cmap=plt.cm.Greys)
ax.locator_params(nbins=3)
ax.set_yticks(np.arange(len(MWQ_keys)))
ax.set_yticklabels(MWQ_keys, fontsize='x-large')
ax.set_xticks(np.arange(3))
ax.set_xticklabels(['RTS', 'BOOTS', 'Common', ], rotation=90, fontsize='large')
ax.vlines(1.5, -0.5, 25.5)
ax.vlines(0.5, -0.5, 25.5)
ax.hlines(12.5, -0.5, 2.5)
divider = make_axes_locatable(plt.gca())
cax_a = divider.append_axes("right", "20%", pad="20%")
cb_a = plt.colorbar(im_a, cax=cax_a)
cb_a.set_label('Canonical weights')
# plt.tight_layout()
In [7]:
matplotlib.rcParams.update({'font.size': 14})
multi_com1 = np.sqrt(corr_mat_1_flat[:, 16 : 32] * corr_mat_1_flat[:, 32:])
multi_com2 = np.sqrt(corr_mat_2_flat[:, 16 : 32] * corr_mat_2_flat[:, 32:])
bin_comp1 = np.array(multi_com1 > 0.1, dtype=np.int)
bin_comp2 = np.array(multi_com2 > 0.1, dtype=np.int)
bin_comp_brain = np.vstack((bin_comp1, bin_comp2))
multi_com1 = np.sqrt(y_loadings_1[:, 1] * y_loadings_1[:, 2])
multi_com2 = np.sqrt(y_loadings_2[:, 1] * y_loadings_2[:, 2])
bin_comp1 = np.array(multi_com1 > 0.1, dtype=np.int)
bin_comp2 = np.array(multi_com2 > 0.1, dtype=np.int)
bin_comp_thought = np.hstack((bin_comp1, bin_comp2))
brain_mat = np.hstack((x_loadings[:, 16:], bin_comp_brain))
behav_arr = np.hstack((y_loadings[:, 1:], bin_comp_thought.reshape((26,1))))
fig = plt.figure(figsize=(17,10))
ax1 = plt.subplot2grid((4,4), (0, 0), colspan=3, rowspan=3)
ax2 = plt.subplot2grid((4,4), (0, 3), colspan=1, rowspan=3)
RSC_plot(brain_mat, ax1)
MWQ_plot(behav_arr, ax2)
for i in range(3):
ax1.annotate('PMC', xy=(0.045 + i*.335, 1.07), xytext=(0.045 + i*.335, 1.1), xycoords='axes fraction',
fontsize='x-large', ha='center', va='bottom',
arrowprops=dict(arrowstyle='-[, widthB=1.2, lengthB=0.6', lw=1.0, color='royalblue'))
ax1.annotate('TPJ', xy=(0.125 + i*.335, 1.07), xytext=(0.125 + i*.335, 1.1), xycoords='axes fraction',
fontsize='x-large', ha='center', va='bottom',
arrowprops=dict(arrowstyle='-[, widthB=1.2, lengthB=0.6', lw=1.0, color='royalblue'))
ax1.annotate('dm\nPFC', xy=(0.185 + i*.335, 1.07), xytext=(0.185 + i*.335, 1.1), xycoords='axes fraction',
fontsize='x-large', ha='center', va='bottom',
arrowprops=dict(arrowstyle='-[, widthB=.4, lengthB=0.6', lw=1.0, color='royalblue'))
ax1.annotate('vm\nPFC', xy=(0.27 + i*.335, 1.07), xytext=(0.27 + i*.335, 1.1), xycoords='axes fraction',
fontsize='x-large', ha='center', va='bottom',
arrowprops=dict(arrowstyle='-[, widthB=2.0, lengthB=0.6', lw=1.0, color='royalblue'))
for i in range(2):
ax1.annotate('PMC', xy=(-0.04, 0.94 - i*.5), xytext=(-0.08, 0.917 - i*.5), xycoords='axes fraction',
fontsize='x-large', ha='right', va='bottom',
arrowprops=dict(arrowstyle='-[, widthB=1.2, lengthB=0.6', lw=1.0, color='royalblue'))
ax1.annotate('TPJ', xy=(-0.04, 0.82 - i*.5), xytext=(-0.08, 0.797 - i*.5), xycoords='axes fraction',
fontsize='x-large', ha='right', va='bottom',
arrowprops=dict(arrowstyle='-[, widthB=1.2, lengthB=0.6', lw=1.0, color='royalblue'))
ax1.annotate('dmPFC', xy=(-0.04, 0.72 - i*.5), xytext=(-0.08, 0.697 - i*.5), xycoords='axes fraction',
fontsize='x-large', ha='right', va='bottom',
arrowprops=dict(arrowstyle='-[, widthB=.4, lengthB=0.6', lw=1.0, color='royalblue'))
ax1.annotate('vmPFC', xy=(-0.04, 0.59 - i*.5), xytext=(-0.08, 0.567 - i*.5), xycoords='axes fraction',
fontsize='x-large', ha='right', va='bottom',
arrowprops=dict(arrowstyle='-[, widthB=1.8, lengthB=0.6', lw=1.0, color='royalblue'))
compNames = ['Positive\nHabitual', 'Spontaneous\nOff-Task']
ax1.annotate(compNames[i], xy=(-0.25, (0.7 - i*0.6)), xytext=(-0.25, (0.7 - i*0.6)), xycoords='axes fraction',
fontsize='x-large', ha='center', va='bottom', rotation=90)
for i, lab in enumerate(['RTS', 'BOOTS', 'Common', ]):
ax1.annotate(lab, xy=((0.5 + i)/3, 1.25), xytext=((0.5 + i)/3, 1.25), xycoords='axes fraction',
fontsize='xx-large', ha='center', va='bottom')
# plt.tight_layout()
plt.savefig('Results/Fig3c_SCCA_components_.png', dpi=300, transparent=True, bbox_inches='tight')
plt.show()
In [379]:
# save for brainview
np.savetxt('comp1_brain.txt', corr_mat_1[..., 2], delimiter='\t', fmt='%10.5f')
np.savetxt('comp2_brain.txt', corr_mat_2[..., 2], delimiter='\t', fmt='%10.5f')
In [393]:
# thought probe
for i in range(2):
behav_arr = np.zeros((13, 1))
behav_arr.flat = y_loadings[0 + 13 * i : 13 * (i +1) , -1]
fig = plt.figure()
ax = fig.add_subplot(111)
im = ax.matshow(behav_arr, vmin=-0.9, vmax=0.9, interpolation='nearest', cmap=plt.cm.RdBu_r, aspect='equal')
divider = make_axes_locatable(plt.gca())
cax = divider.append_axes("right", "30%", pad="20%")
plt.colorbar(im, cax=cax)
ax.set_yticks(np.arange(len(keys)))
ax.set_yticklabels(keys, fontsize='x-large')
ax.set_xticks(np.arange(1))
ax.set_xticklabels('')
plt.savefig('Results/CCA_thought_%i.png' %(i+1), dpi=300, transparent=True, bbox_inches='tight')
plt.show()
In [4]:
region_labels_fn = 'data/data_cross_corr_Bzdok_DMN16_preprocessed_ROIS.npy'
behavefn = 'data/data_MWQ_session_preprocessed.npy'
rscorrfn = 'data/data_cross_corr_Bzdok_DMN16_preprocessed.npy'
import numpy as np
import matplotlib.pyplot as plt
import matplotlib
from os.path import expanduser
matplotlib.rcParams.update({'font.size': 14})
subject_subset = np.load(expanduser(behavefn))[:, 0].astype('i4') - 1
rest_data = np.load(expanduser(rscorrfn))
X = rest_data[subject_subset,:]
mean_NetworkStrength = np.mean(X, axis=0)
n_areas= 16
idx = np.triu_indices(n_areas, 1)
corr_mat = np.zeros((n_areas, n_areas))
corr_mat[idx] = mean_NetworkStrength
corr_mat = corr_mat + corr_mat.T
region_labels = np.load(region_labels_fn)
plt.matshow(corr_mat, vmin=-0.7, vmax=0.7, cmap=plt.cm.RdBu_r)
plt.plot([-0.5, n_areas-0.5], [-0.5, n_areas-0.5], ls='--', c='.3')
plt.xticks(range(16), region_labels, rotation=90)
plt.yticks(range(16), region_labels)
plt.colorbar()
plt.savefig('./Results/RemoveHighFD/average_NS.png', bbox_inches='tight', dpi=150)
In [8]:
import pandas as pd
In [69]:
master_df = pd.read_csv('./data/revision2.csv', na_values=' ')
In [70]:
all_labels = list(master_df.columns)
all_values = np.float64(np.array(master_df.values.tolist()))
var_corr_mat = np.corrcoef(all_values.T)
In [71]:
matplotlib.rcParams.update({'font.size': 14})
plt.matshow(var_corr_mat, vmin=-1, vmax=1, cmap=plt.cm.RdBu_r)
# plt.plot([-0.5, n_areas-0.5], [-0.5, n_areas-0.5], ls='--', c='.3')
plt.xticks(range(var_corr_mat.shape[0]), all_labels, rotation=90)
plt.yticks(range(var_corr_mat.shape[1]), all_labels)
plt.colorbar()
plt.savefig('./Results/RemoveHighFD/wth_2.png', bbox_inches='tight', dpi=150)